Fitting
###Parameter Analyses###
rm(list=ls())
library(ggplot2)
library(mixRSVP)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(magrittr)
library(BayesFactor)
## Loading required package: coda
## Loading required package: Matrix
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************
participantPlots = TRUE #plot histograms with density?
inclusionBF <- function(priorProbs, variable){
###https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp###
if(typeof(priorProbs) == 'S4') priorProbs <- as.vector(priorProbs)
theseNames <- names(priorProbs)
nProbs <- 1:length(priorProbs)
variableMatches <- grep(variable, theseNames)
if(grepl(':', variable)){
subordinateVariables <- variable %>% strsplit(':') %>% unlist()
thisRegex <- paste0(subordinateVariables,collapse = '.*\\+.*')
subordinateEffects <- grep(thisRegex, theseNames, perl = T)
subordinateEffects <- subordinateEffects[!subordinateEffects %in% variableMatches]
sum(priorProbs[variableMatches])/sum(priorProbs[subordinateEffects])
} else {
interactionMatches <- grep(paste0(variable,'(?=:)|(?<=:)',variable), theseNames, perl = T)
variableMainEffects <- variableMatches[!variableMatches %in% interactionMatches]
otherMainEffects <- nProbs[!nProbs %in% c(variableMainEffects,interactionMatches)]
sum(priorProbs[variableMainEffects])/sum(priorProbs[otherMainEffects])
}
}
allErrors <- read.table('../allErrors.txt', sep='\t', stringsAsFactors = F, header = T)
nReps <- 100
bounds <- parameterBounds()
bounds['precision','upper'] <- 3
runAnyway <- FALSE #If TRUE, fit models regardless of the presence of a parameter file.
plots <- FALSE
nParamFiles <- length(list.files(path = '../../modelOutput/',pattern ='parameterEstimates.*\\.csv',full.names = T)) #How many parameter DFs are saved?
if(nParamFiles>0){
paramFiles <- list.files(path = '../../modelOutput',pattern ='parameterEstimates.*\\.csv',full.names = T) #What are the saved param files?
splits <- strsplit(paramFiles, 'Estimates|(?<=[0-9][0-9])_|\\.csv',perl = T) #split out the date stamps from the param files names
theseDates <- lapply(splits,FUN =function(x) paste(x[2],x[3])) %>%
unlist %>%
as.POSIXct(.,format='%d-%m-%Y %H-%M-%S') #Format them as POSIXct (datetime)
whichMaxDate <- which(theseDates == max(theseDates)) #which is the newest?
params <- read.csv(paramFiles[whichMaxDate]) #load the newest
#Empty space in the param DF for the un-modelled participants
notModelled <- allErrors %>% filter(., !ID %in% params$ID) %>% pull(ID) %>% unique
unModelledRows <- expand.grid(
ID = factor(notModelled),
crowded = factor(unique(allErrors$crowded)),
ring = factor(unique(allErrors$ring)),
efficacy = 999,
latency = 999,
precision = 999,
val = 999,
valGuessing = 999,
pLRtest = 999,
stringsAsFactors = F
)
params <- rbind(params, unModelledRows)
} else {
notModelled <- allErrors %>% pull(ID) %>% unique
params <- expand.grid(
ID = factor(notModelled),
crowded = factor(unique(allErrors$crowded)),
ring = factor(unique(allErrors$ring)),
efficacy = 999,
latency = 999,
precision = 999,
val = 999,
valGuessing = 999,
pLRtest = 999,
stringsAsFactors = F
)
}
if(length(notModelled)>0 & runAnyway){
for(thisParticipant in notModelled){
for(thisCondition in unique(allErrors$crowded)){
for(thisRing in unique(allErrors$ring)){
#print(paste0('Participant: ', thisParticipant, '. Ring: ', thisRing, '. Condition: ', thisCondition))
try({
invisible( #invisible combined with capture output hides those print 'error' statements hidden deep in mixRSVP::
capture.output(
theseParams <- allErrors %>% filter(., crowded == thisCondition, ring == thisRing, ID == thisParticipant) %>% analyzeOneCondition(., 24, bounds, nReps)
)
)
}, silent = T) #Doing everything I can to silence those messages
if(theseParams$pLRtest<.05){
params %<>%
mutate(efficacy=replace(efficacy, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$efficacy)) %>%
as.data.frame()
params %<>%
mutate(latency=replace(latency, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$latency)) %>%
as.data.frame()
params %<>%
mutate(precision=replace(precision, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$precision)) %>%
as.data.frame()
} else {
params %<>%
mutate(efficacy=replace(efficacy, ID == thisParticipant & crowded == thisCondition & ring == thisRing, 0)) %>%
as.data.frame()
params %<>%
mutate(latency=replace(latency, ID == thisParticipant & crowded == thisCondition & ring == thisRing, NaN)) %>%
as.data.frame()
params %<>%
mutate(precision=replace(precision, ID == thisParticipant & crowded == thisCondition & ring == thisRing, NaN)) %>%
as.data.frame()
}
params %<>%
mutate(val=replace(val, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$val)) %>%
as.data.frame()
params %<>%
mutate(valGuessing=replace(valGuessing, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$valGuessing)) %>%
as.data.frame()
params %<>%
mutate(pLRtest=replace(pLRtest, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$pLRtest)) %>%
as.data.frame()
}
}
}
write.csv(params, paste0('modelOutput/parameterEstimates',format(Sys.time(), "%d-%m-%Y_%H-%M-%S"),'.csv'),row.names = F)
}
paramsForAnalysis <- params %>% filter(efficacy>.1 & ID != 'CH')
paramsForAnalysis %<>% filter(!latency >= bounds[2,2] & !latency <= bounds[2,1])
paramsForAnalysis %<>% filter(!precision >= bounds[3,2] & !precision <= bounds[3,1])
paramsForAnalysis$ring %<>% as.factor
paramsForAnalysis %<>% mutate(IDChar=ID) %>% mutate(ID = as.numeric(ID)) %>% mutate(ID=as.factor(ID)) %>% as.data.frame() #For some reason the string 'AC5SONA' is interpreted as a nul by anovaBF, which uses base64
Efficacy Analyses
efficacyBF <- anovaBF(efficacy ~ ring * crowded + ID,
data=paramsForAnalysis,
whichRandom = 'ID'
)
efficacyPriorProbs <- efficacyBF %>% newPriorOdds() %>% `*`(efficacyBF) %>% as.BFprobability() %>% as.vector() #extract P(M|Data) from BF object
efficacyInclusionBFs <- rep(-999, times = 3)
names(efficacyInclusionBFs) <- c('crowded','ring','crowded:ring')
for(name in names(efficacyInclusionBFs)){ #calculate inclusion BFs for main effects and interactions
thisInclusionBF <- inclusionBF(efficacyPriorProbs, name)
efficacyInclusionBFs[name] <- thisInclusionBF
}
knitr::kable(efficacyInclusionBFs)
| crowded |
0.4620900 |
| ring |
135.8714676 |
| crowded:ring |
0.2057757 |
#Only evidence for an effect of ring
ggplot(paramsForAnalysis, aes(x=crowded, y = efficacy))+
geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
geom_point(aes(group = factor(ring)), position = position_dodge(.9))+
stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = efficacy))+
geom_point(aes(colour = ID), position = position_dodge())+
geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
#stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
#stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
lims(y=c(0,1))+
facet_wrap(~crowded)
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Width not defined. Set with `position_dodge(width = ?)`

Latency analysis
latencyBF <- anovaBF(latency ~ ring * crowded + ID,
data=paramsForAnalysis,
whichRandom = 'ID'
)
latencyPriorProbs <- latencyBF %>% newPriorOdds() %>% `*`(latencyBF) %>% as.BFprobability() %>% as.vector() #that nested newPriorOdds() function call is the shittiest hack. Prior odds are 1 by default. I shouldn't need it, but S4 is a mystery to me
latencyInclusionBFs <- rep(-999, times = 3)
names(latencyInclusionBFs) <- c('crowded','ring','crowded:ring')
for(name in names(latencyInclusionBFs)){
thisInclusionBF <- inclusionBF(latencyPriorProbs, name)
latencyInclusionBFs[name] <- thisInclusionBF
}
knitr::kable(latencyInclusionBFs)
| crowded |
2.2341273 |
| ring |
0.4563018 |
| crowded:ring |
0.5005254 |
ggplot(paramsForAnalysis, aes(x=crowded, y = latency))+
geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
geom_jitter(aes(group = factor(ring)), position = position_dodge(.9))+
stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = latency))+
geom_point(aes(colour = ID), position = position_dodge())+
geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
#stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
#stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
facet_wrap(~crowded)
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Width not defined. Set with `position_dodge(width = ?)`

Precision Analysis
precisionBF <- anovaBF(precision ~ ring * crowded + ID,
data=paramsForAnalysis,
whichRandom = 'ID'
)
precisionPriorProbs <- precisionBF %>% newPriorOdds() %>% `*`(precisionBF) %>% as.BFprobability() %>% as.vector() #that nested newPriorOdds() function call is the shittiest hack. Prior odds are 1 by default. I shouldn't need it, but S4 is a mystery to me
precisionInclusionBFs <- rep(-999, times = 3)
names(precisionInclusionBFs) <- c('crowded','ring','crowded:ring')
for(name in names(precisionInclusionBFs)){
thisInclusionBF <- inclusionBF(precisionPriorProbs, name)
precisionInclusionBFs[name] <- thisInclusionBF
}
knitr::kable(precisionInclusionBFs)
| crowded |
0.3971054 |
| ring |
2192.4479671 |
| crowded:ring |
0.2211942 |
ggplot(paramsForAnalysis, aes(x=crowded, y = precision))+
geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
geom_jitter(aes(group = factor(ring)), position = position_dodge(.9))+
stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = precision))+
geom_point(aes(colour = ID), position = position_dodge())+
geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
#stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
#stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
facet_wrap(~crowded)
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Width not defined. Set with `position_dodge(width = ?)`

Plots with Density
paramsForAnalysis %<>% mutate(ID = IDChar) %>% as.data.frame #Convert the problem name back for plotting
if(participantPlots){
for(thisParticipant in unique(paramsForAnalysis$ID)){
theseErrors <- allErrors %>% filter(ID==thisParticipant)
theseDensities <- data.frame(density = numeric(0),SPE = numeric(0), crowded = character(0), ring = numeric(0), stringsAsFactors = F)
theseConditions <- paramsForAnalysis %>% filter(ID == thisParticipant) %>% pull(crowded) %>% unique
for(thisCondition in theseConditions){
theseRings <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition) %>% pull(ring) %>% unique
for(thisRing in theseRings){
thisEfficacy <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(efficacy)
thisLatency <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(latency)
thisPrecision <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(precision)
thisConditionErrors <- theseErrors %>% filter(crowded == thisCondition & ring == thisRing)
# print(paste0('Participant: ', thisParticipant, '. Ring: ', thisRing, '. Condition: ', thisCondition,'. N = ', nrow(theseErrors)))
minError <- thisConditionErrors %>% pull(SPE) %>% min
maxError <- thisConditionErrors %>% pull(SPE) %>% max
thisRange <- seq(minError,maxError,.1)
thisConditionDensities <- data.frame(SPE = thisRange, density = dnorm(thisRange, thisLatency, thisPrecision), crowded = thisCondition, ring = thisRing, stringsAsFactors = F)
theseDensities <- rbind(theseDensities, thisConditionDensities)
# if(any(is.nan(theseDensities$density))){
# print(theseDensities)
# }
# thisFileName <- paste0('modelOutput/Plots/',thisCondition,'/',thisRing,'/',thisParticipant,'-',format(Sys.time(), "%d-%m-%Y_%H-%M-%S"),'.png')
# ggsave(filename = thisFileName, thisPlot,width = 16, height = 9)
}
}
thisPlot <- ggplot(theseErrors, aes(x=SPE))+
geom_histogram(binwidth = 1)+
geom_line(data = theseDensities, aes(x = SPE, y=density*60))+ #scale density to histogram with density * N * binwidth
scale_y_continuous(sec.axis = sec_axis(~./60, name = 'Density'))+
labs(y = 'Frequency', title = thisParticipant)+
facet_wrap(~crowded*ring)
show(thisPlot)
}
}














